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Abstract. We study, experimentally and theoretically, the controlled transfer of 
harmonically trapped ultracold gases between different quantum states. In particular 
we experimentally demonstrate a fast decompression and displacement of both a non- 
interacting gas and an interacting Bose-Einstein condensate which are initially at 
equilibrium. The decompression parameters are engineered such that the final state is 
identical to that obtained after a perfectly adiabatic transformation despite the fact 
that the fast decompression is performed in the strongly non-adiabatic regime. During 
the transfer the atomic sample goes through strongly out-of-equilibrium states while 
the external confinement is modified until the system reaches the desired stationary 
state. The scheme is theoretically based on the invariants of motion and scaling 
equations techniques and can be generalized to decompression trajectories including 
an arbitrary deformation of the trap. It is also directly applicable to arbitrary initial 
non-equilibrium states. 



PACS numbers: 67.85.-d, 37.10.-x 
Introduction 

In Quantum Mechanics, the evolution of a system described by a time-dependent 
Hamiltonian H(t) is adiabatic when the transition probabilities between the 
instantaneous eigenstates of H are negligible. This happens when H is either time- 
independent, or when its rate of change is slow compared to the typical time-scales 
involved [IH3]. Nevertheless, thinking in terms of instantaneous eigenstates is often 
much easier than looking for the solutions of time-dependent problems. In the field 
of atomic physics, going from the semi-classical approach of atom-field interaction to 
the celebrated dressed state picture [I] illustrates the convenience of such adiabatic 
representations. 

For this reason, many adiabatic schemes to prepare interesting quantum states 
were proposed. For instance, non-classical states [HI E], or new strongly correlated 
states [7] can be prepared by adiabatic passage. Quantum adiabatic computation 
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has recently been demonstrated [8]. Yet adiabatic techniques are typically slow [3], 
while experimentalist are often constrained by finite lifetimes or coherence times of 
their samples. This motivated the search for fast schemes reproducing or approaching 
adiabatic transformations. Some methods use minimization techniques to optimize the 
transition to a target state [9HT2] , whereas others yield the exact same state that would 
have been reached after an adiabatic transformation [131 HI]- The latter are referred 
to as shortcuts to adiabaticity. In this article, we detail how such methods can be 
used on the motional degrees of freedom of ultracold gases confined in time-dependent 
harmonic traps, and experimentally demonstrate the validity of the approach. Two 
direct applications of the procedure are the fast cooling of atomic samples, and the 
suppression (or reduction) of any parasitic excitations which occur in experiments on 
ultracold gases when the trap geometry or the interactions are modified. Since the 
method is not restricted to equilibrium states it could be used in a variety of situations 
as discussed at the end of the paper. 

The first part is theoretical and recalls how harmonically confined gases react 
to the variation of the trap. Both the one- dimensional non-interacting gas, and the 
three-dimensional Bose-Einstein condensate with repulsive contact interaction between 
particles are treated. In the second part, the method to realize shortcuts to adiabaticity 
are detailed for theses two systems, and examples are given. The third part focuses 
on the experimental realization of these methods. Rapid decompressions have been 
performed on both a non-interacting gas and a Bose-Einstein condensate. The practical 
limitations which degrade the results are discussed. In the last part of the article, we 
attempt to generalize the problem to an arbitrary variation of the three-dimensional 
harmonic potential and give other examples of shortcuts which may be of experimental 
relevance. 

1. Scaling properties of harmonically confined ultracold gases: two 
examples 

In this section, we recall how the density and velocity distributions of a one-dimensional 
(ID) non-interacting gas are affected by a change of the harmonic confinement. In ID, 
the harmonic trap is fully described by its time-dependent angular frequency u(t), and 
minimum position qo(t). We show that the dynamics is fully described by two scaling 
functions, one associated to the cloud's size, the other to its centre-of-mass position, and 
exhibit the exact solutions of the Schrodinger equation. This will be used in the rest 
of the paper to realize shortcuts to adiabaticity (cf. Sec. |2|. Similar scaling properties 
are also recalled for Bose-Einstein condensates (BECs) with strong interactions in the 
Thomas-Fermi regime. The analogy between the invariant method used for the non- 
interacting gas [To] , and the scaling often used for BECs [T6T4T8"] is underlined. 
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1.1. The non-interacting gas 

We consider a ID non-interacting gas confined in the most general time-dependent 
harmonic potential, described by the one particle Hamiltonian 

v 2 1 

H(q,p,t) = j- + -mu 2 (t) [q - q (t)] 2 , (1) 

where q and p are conjugate variables, and m is the mass of a particle. We first recall 
how dynamical invariants can be used to find the general solutions of the Schrodinger 
equation. 



1.1.1. Definition and properties of dynamical invariants In 1969 Lewis and Riesenfeld 
[15] generalized the concept of invariant of motion to the case of explicitly time- 
dependent Hamiltonians H(q,p,t). Such Lewis invariants (also called dynamical 
invariants, or first integrals) can be used to solve the Schrodinger equation 

ih^ = H{q,p,t)\t). (2) 

Given a time-dependent Hamiltonian H(q,p,t), a time- dependent hermitian 
operator I(q,p, t) is a dynamical invariant of the system described by H if it is constant 
under Hamiltonian evolution, that is if 

dt dt ih> ' J K ' 

In this case, the following properties hold [15J: 

(i) if \t) is a solution of (j2l), then I\t) is also a solution of ([2j, 

(ii) the eigenvalues A(t) and associated eigenvectors \X;t) of / are a priori time- 
dependent. We assume they form a complete set. It turns out that the eigenvalues 
are actually constant: \(t) = A. They are real because I is hermitian. 

(iii) The eigenvectors of I satisfy 

d 

for all A, A' such that A ^ A', (A'; t\ih~ — \\; t) = (A'; t\H\\; t). (4) 

(iv) If we assume that / does not contain the operator d/dt, then the set of vectors 
{e tax ^\\; t), a\(t) G is also a complete set of eigenvectors of /. If these 
functions are chosen to solve the equations 

dct\ , . d H . , , 

■a-HM-s-alA;*) («) 

then Eq. Q4J) also holds for A' = A. Using the fact that the set is complete, this 
gives the general solutions of the time-dependent Schrodinger equation as 

|t) = ^c A e^W|A;t), (6) 

A 

where the ca's are constant complex numbers. 
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The solutions of the Schrodinger equation are thus given by the knowledge of an invariant 
I(q,p,t), any set of its time-dependent eigenvectors, and the phases a\(t) which must 
solve Eqs. 

1.1.2. Derivation of a dynamical invariant In this section, we give a simple derivation 
of the invariants of a ID time- dependent harmonic oscillator (HO) described by Q. 
We use the classical formalism to derive the invariant, which is also an invariant of the 
corresponding quantum system. 

The canonical equations of motion associated with the Hamiltonian ([I]) are 

^ = {p,H} = -mu> 2 (t)[q-q (t)], (7b) 

i r a m dAdB dAdB . 

where {A, i?} = — — — — are the roisson brackets of two observables A and B. 

dq op op dq 

When the angular frequency u(t) and trap centre qo(t) vary, one expects the cloud 
to be displaced and to change its size, thus one can introduce a canonical change of 
variables 

Q= q ~ b lf ] > P = P(Q,P,tl r = r(t), (8) 

leading to a new Hamiltonian H'. One has to derive conditions on the real dimensionless 
function b, and the displacement function q cm such that the transformation is canonical. 
For this, we look for a new Hamiltonian of the form 

H' = ^ + \mulQ^ + f{r), (9) 

where Uq is a constant angular frequency. The Hamiltonian explicitly depends on time 
only through the function /(r) (which does not contain the variables Q and P). The 
transformation ^ is canonical if 

d Q = {Q,H>), (10a) 



From Eq. (10a) one deduces that 
and that 



dr = b- 2 dt (11) 
P = b(p-mq cm ) -mb(q- g cm ), (12) 



where ' denotes the derivation with respect to time t. From Eq. (10b), one finds the 
functions b and g cm must obey the two differential equations 

b + uj 2 (t)b=^, (13) 

q cm + co 2 (t)[q cm (t) - q (t)} = 0. (14) 
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When these two equations are satisfied, the quantity 

/= £ + ? n ^ 2 (i5 » 

which appears in the new Hamiltonian is a Lewis invariant. This can be proved directly 
by checking that Eq. ^ is verified. 

The choice of the function /(r) in H' is irrelevant for the dynamics, since doing the 
change of Hamiltonian 

H'^H'- f(r) = I (16) 

corresponds to a gauge transformation which changes the phase of the wave function in 
the following manner: 

1>H>(Q,T)^MQ,T) = e* F V>t/> H ,{Q,T), (17) 
where F is a primitive of /. 



1.1.3. Wave functions Once an invariant has been found, the results of section 1.1.1 
can be used to calculate the wave functions of the time-dependent HO ([!]). We use 
Dirac's method to calculate the wave function of the time-independent HO ( Jl5j ). We 
define dimensionless variables 



V n ymnujQ 

satisfying the commutation relation [£, it] = i, and introduce the lowering and raising 
operators 

a=-^(£ + i7c), a* = -j=(Z-m). (19) 

The invariant reads 

I = hu {a i a + 1/2). (20) 
The eigenstates |n) of the number operator h = a) a are the eigenstates of I and satisfy 

a\n) — y/n\n — 1), a)\n) = y/n + l\n + 1), nGN. (21) 

The eigenvalues of I are 

A„ = (n + M hu , neN. (22) 
The wave function il> (q,t) = (q\0) is calculated by solving the equation 

o|0) = (23) 



in \q) representation. The expression of it is obtained from p = —ihd/dq, and Eqs. (18) 
and (J8J), and reads 

d bb _ / m , . . 

tt = -z— i-\h—bq cm . 24 
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Imposing the normalization condition J dq \ipo{q, t)\ 2 = 1, and calculating the time- 
dependent phase corresponding to the initial Hamiltonian 0, we obtain the wave 
function 



7T 



-1/4 



exp 



1 ( q - q CI 

2 V ahob 



-£f(t) Mq,t) p-iAo/o* dt'/ft 3 



(25) 



where 



</>(q,t) 



m 



+ ^ ( <knJ> - <In„h j <l 



F(t) 



m 
2~ 



dt' 



2y C m , ,2/./n 2 
W 0"^T + 6(7 l J J% 



(26) 
(27) 



and g 0) Q'cm, and their derivatives are functions of t (t' when they are under an integral 
symbol) and are linked by Eqs. (13) and (14). a^ = ^h/mujQ is the HO length of /. 

From this expression, we see the physical meaning of the two scaling functions: 
9cm (t) is the centre of the wave function (centre of mass of a cloud which was initially 
at equilibrium), and a^b is the width of the wave function. 

The wave function associated to the eigenvalue A n of I is expressed in terms of the 
nth Hermite polynomial H n as 



ipn(q,t) 



2 n n\ 



i/j (q,t)H n 



q-q C r 

a ho b 



-i(\n-\ )tidt'/b 2 



(2c 



1.2. The case of an interacting Bose-Einstein condensate 

For the corresponding three-dimensional (3D) interacting system of iV particles, the 
Hamiltonian is 

,2 



H 



N 

E 

i=i L 



2m 



+ U{r t ,t) 



E 

i<3 



V(rj - ri 



(29) 



The potential U is supposed to be a time-dependent 3D HO, and the rotation of this 
harmonic confinement is excluded for the moment (the trap keeps the same eigenaxes): 



U(r,t) = -mU(t) [r, - rl(t)] 2 + «J(f) [r y - r° y (t)] 2 + u 2 z (t) [r z -rl{t)\ 



(30) 



V is the interaction potential between two particles, which is well approximated by a 
delta function for ultracold gases [19|. 



The procedure described in Sec. |1.1| cannot be easily adapted, because it would 
require the knowledge of an invariant of this many-body system. But, when dealing 
with a BEC, the dynamics is well described by a single particles wave function, 
whose evolution obeys a non-linear Schrodinger equation, the Gross-Pitaevskii equation 
(GPE) 
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1.2.1. Scaling approach Let us consider a quantum system described by the wave 
function ip(r,t), whose time evolution is given by the GPE 



_|l A + f/M) + WV|^(r,t)| 2 
2m 



^(r,t), 



(31) 



with m, the mass of the particles, N the number of particles, and V = Aith 2 a s /m 
the interaction coupling constant generated by s-wave scattering between particles, 
characterized by the scattering length a s . Analogously to the non-interacting case, 
we wish to write the solution of the time-dependent GPE as a function of the solution 
of a time- independent one expressed in a suitable frame of reference. Following this line, 



a strategy to solve Eq. (31) is to find a change of variables p(r, {bi(t)}, {^ m (t)}) where 
the biS and the rf m 's are scaling and translation functions such that Eq. (31) can be 



written as a time-independent equation (i.e. a GPE with a time-independent potential) 
on the wave function x(Pj t )> defined by the relation 



V>M) =A(t) X (p,r)e l 



(r,t) 



(32) 



A(t) being a time-dependent normalization factor and (j)(r,t) a space- and time- 
dependent phase. All the dynamics induced by the time-dependent potential is 
transferred to the functions {bi(t)} and {r^ m (t)}, and the differential equations they have 



to satisfy (similar to Eqs. (13) and (14)). If one can solve the new time-independent 



equation on x> one solves Eq. (31) and knows the wave function ij)(r,t). 



Equation (31) is invariant under the transformation 

\/ie{x,y,z}, pi = 



(t) 



bi(t) 



(33) 



in any of the following cases: 
(i) in the non-interacting limit [16, 20J: in this case the system is equivalent to three 



independent HO of the kind treated in Sec. 1.1 



(ii) for a suitable driving of the interaction term V [20], that is, assuming one can 
control V(t) at will (for cold gases, this can be done using Feshbach resonances), 

(iii) in the TF limit [IT]. 

This third case, which is detailed in the following section, is used in the rest of the 
paper. 



1.2.2. Condensate wave function in the Thomas-Fermi approximation Given a time- 
dependent Gross-Pitaevskii equation, the TF approximation consists in neglecting the 
kinetic-energy-like term in the p-frame of reference, i.e. —h 2 /{2m) b~[ 2 d 2 x/ dp 2 , 
supposed to be small compared to the interaction term [TTJEIE]- In this regime, provided 
that A{t) = (ILA)~ 1/2 and that 



M) 




2 h h 



(34) 
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with 

Mi) 



rn 

2h 



cm\ 2 



where the scaling and translation functions are solutions of 

^ 2 (o) 



Vz G bi + lo { (t)bi 



kb x byb z ■ 



rT + ^it) [rT "»?(*)] =0, 
one gets the following equation on 

ih^ X {p,r) = [u(p,0) + VN\x(p,t)\ 2 ] x (p,r) 



where we defined a rescaled time 



r(t) 



n f 6i(f) 



If at t — the condensate is at equilibrium, the solution of Eq. (38) is 



VN 



(36) 
(37) 

(38) 

(39) 



(40) 



/j, being the chemical potential. This gives the typical inverted parabola density profile 
whose sizes evolve in time as Ri(t) = Ri(0)bi(t), Ri(0) = ^/2jj, /mujf(0) being the initial 
TF radii. 



2. Shortcuts to adiabaticity 

In this section the definition of a shortcut to adiabaticity is given, and the results of 
Sec. [T] are used to derive angular frequency trajectories realizing such shortcuts, for both 
non-interacting gases and interacting BECs confined in time-dependent harmonic traps. 



2.1. Shortcut to adiabaticity based on an invariant of motion 

For a system described by a Hamiltonian H(t), a shortcut to adiabaticity is realized 
when another Hamiltonian H'(t) can be found, such that the state obtained after a 
finite time of evolution with H'(t) is identical (up to a global phase factor) to the final 
state of the adiabatic evolution with H(t). Shortcuts to adiabaticity are not adiabatic; 
only the final state is identical to that obtained after an adiabatic evolution. 

The possibility of such shortcuts has been known for a long time. For instance, 
in the case of a HO with a time-dependent frequency u(t) treated in Ref. |15j . when 
discussing the transition probability P sm between two instantaneous eigenstates \s;t) 
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and \m ; t), the authors noticed that some trajectories u(t) could lead to the same result 
as the adiabatic case, namely 

Psm 0~sm- (^1) 

Such shortcuts to adiabaticity can thus be realized simply by engineering the time- 
dependent parameters of the Hamiltonian. 

A practical method to find a class of appropriate u(t) was detailed by Chen et al. 
[14J. In this case, the Hamiltonian is chosen to be time-independent (but with different 
frequencies) outside the time interval t G [0, t/]. An invariant is engineered to commute 
with the Hamiltonian outside this interval. This yields a specific u(t) for which all the 
eigenstates of H(t < 0) are exactly mapped to the corresponding ones of H(t > tf) 
after the evolution for t G [0, £/]. Up to a global phase and a rescaling of energies and 
lengths, the final state (at time t = tf) is identical to the initial one (t = 0), i.e. if the 
initial state was 

\ifj ■ t < 0) = °n\ n I t = 0)e-^" (0) *, (42) 

n 

where {\n; t),n 6 N} is a basis of instantaneous eigenstates of H(t) and {hoj n (t)} the 
corresponding eigenvalues, and J^ ra |c n | 2 = 1, then the final state is 

\ip; t > t f ) = e i$ ^c n |n; t f )e~^ {t ^. (43) 

n 

This is true even if the initial state is not an equilibrium state. 



2.1.1. Frequency trajectory for a non-interacting gas The Hamiltonian is assumed to 
have the form 

p 2 1 

H = — + -m tu 2 (t)q 2 + mgq, (44) 

which is identical to Q, with the additional constraint qo(t) = —g/uo 2 (t) (and a gauge 
transformation consisting in adding — m u {t) ql(t) / 2 to H). It describes a single particle 
in a harmonic trap subject to a constant force, which, in the experiments presented in 
Sec.|3j comes from gravity. The angular frequency u(t) is assumed to be constant outside 
the interval t G [0, tf]. During this interval, the problem is to find the appropriate 
frequency trajectory u(t), connecting the initial trap of initial frequency oo(0) to a final 
trap of frequency co(tf), for the decompression (or compression if u(0) < uj(tf)) to 
implement a shortcut to adiabaticity. Figure [l] shows the initial and final situations 
assuming a decompression (co(tf) < u>(0)). 

We used the strategy introduced by Chen et al. [H]. If the invariant commutes 
with the Hamiltonian 

[I,H\ = (45) 

for t < and t >tf, and provided that the functions b and g cm are sufficiently continuous, 
the stationary states of H(t < 0) will be transferred to the corresponding ones of 
H(t > t f ). 



Shortcuts to adiabaticity for trapped ultracold gases 



10 



magnetic 

V potential U 




/ 



co(0) (o(t f ) 



Figure 1. Schematic representation of the trap decompression. The potential (plain 
blue line) is the sum of the gravitational potential (dashed black line) and the harmonic 
magnetic potential (dashed red line). When the trap frequency is changed from u>(0) to 
u>(tf), the lengths are multiplied by 7 = y/oj(0)/u}(tf), and the energies divided by 7 2 . 
Because of gravity, the trap centre shifts vertically by Aq — —g \l/u) 2 (tf) — l/u; 2 (0)] . 



It is convenient to use the dimensionless function 

c W = -f^f (46) 



instead of q cm , and to rewrite Eq. ( 14 ) using the rescaled time r (Eq. (jTTj) ) . Equation ( 14 ) 
becomes 

d 2 c/dr 2 + ul c = tu 2 b 3 . (47) 
If one chooses to set uo = w(0), and the conditions 

6(0) = 1, 6(0) = 0, (48a) 

c(0) = 1, 6(0) = 0, (48b) 

then 1(0) = H(t < 0), and if 



b(t f ) = 7, b(t f ) = 0, (48c) 

c(t f ) = 7 3 , c(t f ) = 0, (48d) 

where 7 = \Ju$/uf, then I (if) = 7 2 i?(t > tf) +h(t), where h is a function of time only. 



These boundary conditions thus fulfil the condition (45). Since the functions b and c 



must be solutions of Eqs. (13) and (47), four additional boundary conditions must be 
satisfied: 

6(0) = 0, b(t f ) = 0, (48e) 

c(0) = 0, c(t f ) = 0. (48f) 
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In order to construct the functions b and c satisfying these boundary conditions and 



the two differential equations (13) and (47), it is convenient to write all the boundary 



conditions on the function c and its derivatives with respect to the rescaled time r. 



Using Eqs. (11) and (13), and differentiating Eq. (47) twice with respect to r, one finds 



the ten conditions 



c(0) = 1, 
c(r/) = 7 3 , 



(49a) 
(49b) 



and, for all k G {1,2,3,4}, 



d k c 
dr k 
d k c 



(0) = 0, 



dr k 



(Tf) = 0, 



(49c) 
(49d) 



which are sufficient for the twelve boundary conditions (48). r/ is the rescaled time 
corresponding to tj: r/ = J^ f b~ 2 (t') dt'. 

Under this form, the boundary conditions are well suited to use a polynomial ansatz 



for c(t), deduce 6(r) with Eq. (47), compute r(t) by numerically integrating Eq. (11), 



and obtain b(t). The final step consists in using Eq. (13) to obtain the time-dependent 
trap frequency as uj 2 (t) = uj\jb^ — bjb. 

An example of this procedure is given on Fig. [2] for particular values of the initial 
and final frequencies. The final rescaled time 7/ can be chosen at will, it can be 
arbitrarily small, but one important constraint on the function c is that it must not 
lead to vanishing values of b which give infinite u 2 (t). Additional constraints on c arise 
from experimental requirements, such as positive oo 2 (t) (attractive potentials), maximal 
and minimal frequencies attainable with a given setup, speed at which u(t) can be varied 
etc. Since all this depends on a particular experimental setup, no mathematical analysis 
of the best ansatz to use was done. 

For the experiments presented in Sec. [3] and in Refs. [2Tj a polynomial of order 
fifteen was used: 

15 

J2 Ck 



c(r) 



k=0 



T 



(50) 



The first coefficient is fixed to 1 by Eq. (49a) and ci, 
We arbitrarily impose C5 



ClO 



• ■ , c 4 are fixed to by Eqs. (49c). 
0, which leaves five coefficients 



which are uniquely determined by the remaining boundary conditions (49b) and (49d). 



The calculation of these remaining coefficient is done by inverting the linear system 
corresponding to these five equations. 

In principle, since there are ten BCs, a 9th order polynomial can be used, which 
yields a unique solution for the ten coefficients of c. Nevertheless, the obtained trajectory 
was not well behaved enough to be realized experimentally (the frequency was decreasing 
too fast in the beginning compared to what could be achieved by the apparatus). This 
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is the reason why a higher order polynomial was used and six coefficients were fixed to 
0. 

Since the polynomial can be of any order greater than 9, and the boundary 
conditions only impose a linear relation between nine of its coefficients, there is obviously 
an infinity of different solutions connecting two given initial and final states. Moreover, 
other functions than polynomials could be used for c, as long as they provide enough 
free parameters. 

The obtained nonzero coefficients of (50) are given in table [I] 



Table 1. Nonzero coefficients of the polynomial ansatz for c(r) calculated from the 



boundary conditions ( 49 ) 



CO 


Cll 


Cl2 


Cl3 


Cu 


Cl5 


1 


1365( 7 3 - 1) 


5005(7 3 - 1) 


6930( 7 3 - 1) 


4290( 7 3 - 1) 


1001(7 3 - 1) 



2.1.2. Example In this section we determine the trajectory used in Sec. |3.2| and in 
Ref. [21] ■ The parameters are given in table [2j Figure [2] shows the functions c(r), b(r), 
t(r) and u)(t)/2ir corresponding to this decompression. 

Table 2. Parameters of the ID decompression of a non-interacting thermal gas. 



Initial frequency uj(0)/2tt 
Final frequency <-o(tf) /2tt 
Final rescaled time Tf 
Corresponding duration tf 



235.8 Hz 
15.67 Hz 
5.65 ms 
35.0 ms 



Since the exact wave functions are known, all the properties of the atomic cloud 
can be calculated during decompression. For instance, Fig. [3] displays the size and 
centre-of-mass position of a cloud initially at equilibrium in the compressed trap. These 
are compared to the same values if the decompression were done very slowly as in the 
adiabatic theorem. The clear difference between the plain and dashed curves illustrates 
the fact that the decompression is not adiabatic. 

2.2. Shortcut to adiabaticity for an interacting Bose-Einstein condensate in the 
Thomas-Fermi limit 



Let us suppose that ip(r,t < 0) is a stationary state of Eq. (31). We can engineer the 



parameters of the potential U (r, t) such that i/j(r,tf) is also a stationary state for t > tf. 



This implies that x(Pi T — r /); with r/ = must be a stationary state of Eq. (38) 

and that V r 0(r,tj) = 0. If these two conditions hold, ip(r,t) can evolve during the 

time interval [0, tf] between two stationary states even being strongly different from the 

adiabatic stationary state during the evolution for < t < tf. In our experiment, the 

time- dependent trapping potential has a cylindrical symmetry of the form 

1 1 
U(r,t) = ^mu 2 ± (t)(r 2 x + r 2 z ) + -mu 2 (t)r 2 y + mgr z , (51) 
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_l i I i I i I i I i I i I U i I i I i I i I i l_ 

012345 012345 
Rescaled time r(ms) r(ms) 




r(ms) £(ms) 



Figure 2. Determination of the frequency trajectory when the trap is decompressed 
from ui(t = Q)/2ir — 235.8 Hz to uj(tf)/2TT = 15.67 Hz within 35 ms (cf. parameters of 
Tab. [2]). (a) A fifteenth order polynomial ansatz is used for the displacement function 
c(t), which gives (b) the scaling function 6(r) through Eq. (47); (c) real time t(r) is 
calculated by numerically integrating Eq. (11); (d) Eq. (13) is used to determine the 
time-dependent frequency w(i)/27r (notice the logarithmic scale). 



with initial and final angular frequencies Co^q^O) and u± t \\(tf) = (0)/7i 11 
respectively. This case corresponds to fix V t, r® y (t) = in Eq. (30) and r° z {t) = 
— g/uj_(t). By introducing the dimensionless function 



c(t) 



"1(0) rT(t) 
9 b ± (t) 



(52) 



the differential equations (36) and (37) take the form 



b ± (t) + b ± (t)ul(t)=ul(0)/[bl(t)b {l (t)] (53) 
6||(t)+6|,(t) W j[(t)=a;|(0)/[6jj(t)^(t)] (54) 
&i (t)&[|(t)c(t) + 2&i (t)&|| (t)b± (t)c(t) + ul (0)c(t) - ul(0)bl(t)b {l (t) = 0. (55) 



The final state is an equilibrium state if the final TF radii verify that 



R±,\\(t f )/R±,\\(0) 
rr(f/)Ar(o) = 



= 7! || , if the vertical centre-of-mass position fulfils the condition 
7Jl, and if the condensate flow is null, namely if V0 = 0. These 
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Figure 3. Expected (a) centre-of-mass position and (b) cloud size during a fast 
decompression (same parameters as Tab. [5] and Fig. [2]). The dashed curves correspond 
to the same values in the adiabatic limit tf — > oo. The adiabatic centre-of-mass 
position is the trap minimum q&d.(t) = — 5/ w2 (^)i an d the adiabatic size is <7ad.(t) = 

y/wo/u)(t)<j{0). 



c{t f ) = b A 

6/5 -2/5 



,11 (0) = Ktf) 



±^ = and c(0) = 1, 

4/5 8/5 



b±(tf) = 7 ± 7|| and b\\(tf) = 7j_* /0 7|| /0 - These latter 
must hold as well, giving sixteen independent boundary 



lead to the boundary conditions c(0) 

imply that b± t \\ (0) = &ui(t/) = 
conditions (BC). 

Our procedure to engineer cu± t \\(t) is to reduce the dimensionality of the problem 
by only considering the trajectories that lead to a constant axial size. This corresponds 



to keeping b\\(t) = b\\(Q) for any t, fixing a trap decompression with j± 
case, Eqs. (|53|-([55| reduce to 

b±(t) 



bi(t)c(t) 



b ± (t)ul(t)=ul(0)/bl(t) 
uj [] (t)=u ll (0)/b ± (t) 
2b 3 ± (t)b ± (t)c(t) + ul(0)c(t) - u> 2 ± (Q)b 3 ± (t) = 0. 



In this 

(56) 
(57) 
(58) 



Equation (56) is identical to Eq. (13) and Eq. (58) is nothing but Eq. (47) expressed 



with the real time (the rescaled time being given by Eq. (39) instead of Eq. (11)). Thus 
we can exploit for b±(t) and c(t) the solutions obtained for the non-interacting gas, 



provided that the axial frequency is varied according to Eq. (57) 



2.2.1. Example As an example of the procedure described above, we determine the 



trajectories used in Sec. 3.3 and in Ref. [22]. The decompression parameters are given 



in table |3j The radial frequency is reduced by a factor of 9, and the axial frequency by 
a factor of 3. 

The obtained trajectories are represented in Fig. |4| 



2.2.2. Validity of the Thomas- Fermi approximation To check the validity of the 
Thomas-Fermi approximation that led to the trajectories of Fig. |4| three-dimensional 
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Table 3. Parameters of the 3D decompression of an interacting Bose-Einstein 

condensate. 

Initial radial frequency wj_(0)/27r 235.8 Hz 

Final radial frequency Lu±(tf)/2ir 26.2 Hz 

Initial axial frequency <jii(Q)/27r 22.2 Hz 

Final axial frequency wii(t/)/27r 7.4 Hz 

Final rescaled time tj 11.555 ms 

Corresponding duration tj 30.0 ms 




DC 



10 20 

Time (ms) 



10 20 
Time (ms) 



Figure 4. (a) Radial and (b) axial trap frequencies for the shortcut decompression of 
a BEC in t f = 30 ms. 



Gross-Pitaevskii simulations have been performed and compared with the analytical 



results of Sec. 2.2 In the numerical solution we use a split step operator in time 
combined with a fast Fourier transformation in space. The results are presented in 
Fig. [5] and show that this approximation is well justified for our experimental parameters 
(decompression of Fig. |4| number of atoms N ~ 10 5 , scattering length of 87 Rb of 
a s ~ 100 a , a being the Bohr radius). 



3. Experimental realization of shortcuts to adiabaticity 

The procedure described above was tested experimentally by quickly decompressing a 
trapped ultracold gas of 87 Rb atoms. In this section, we describe the experimental 
steps involved in the preparation of the cold sample (cold thermal gas or BEC) and 
then explain how the decompression is controlled, monitored and compared to simpler 
(non-optimal) schemes. 

3.1. The apparatus 

The Bose-Einstein condensation apparatus, sketched in Fig. |6j is formed of two chambers 
connected by a differential pumping tube. Each chamber is pumped by a separate ion 
pump. A copper tube containing solid Rubidium provides gaseous Rubidium to the 
upper chamber, resulting in a pressure of ~ 10 -9 mbar (100 nPa) which loads a large 




Figure 5. Comparison between the GPE simulations (dashed red lines) and the 
scaling solutions given by the Thomas-Fermi approximation (solid black lines) showing 
its validity, (a) Centre-of-mass position; (b) axial and radial sizes. The peak relative 
difference between the scaling solution and the GPE simulations are respectively 0.3 % 
and 0.2 % for the axial and radial sizes. The decompression occurs during the first 
30 ms, after which the cloud evolves in the static final trap. 



magneto-optical trap (MOT1). The lower chamber is a glass cell, in which a second 
MOT can be produced. The low-conductance tube connecting the chambers (length 
10 cm, diameter 5 mm) results in a pressure on the order of 10 -11 mbar in the second 
chamber. This low pressure is essential for the production of BECs because background 
gas collisions with the magnetically trapped atoms is the key effect limiting the efficiency 
of evaporative cooling. 

3.1.1. Production of ultracold clouds For the production of a BEC, the first step is the 
loading of MOT2 from MOT1. The light of both traps is red-detuned by 5 = —3.5 T 
from the |5 2 Si /2 , F = 2) |5 2 P 3/2 , F = 3) cooling transition of 87 Rb (r/27r = 6.07 MHz 
is the natural line width of the D2 transition of 87 Rb). The six beams of both traps also 
contain repumper light tuned to the \5 2 Si/2,F — 1} — > \5 2 P^/2,F = 2) transition, which 
prevents atoms from accumulating in the lowest energy state |5 2 ^i /2 ,F = 1). The light 
is provided by two DFB diode lasers, both injected in a single-pass tapered amplifier. For 
both MOTs the light is delivered to the atoms by six polarization-maintaining optical 
fibers, to ensure a good long-term stability of the alignment. The total laser powers in 
MOT1 and 2 are 360 mW and 73 mW respectively, with beam waists of 2.7 cm and 
6.7 mm. The magnetic field gradients of MOT1 and 2 are respectively B[ = 11 G/cm 
and B' 2 = 14.6 G/cm (these values correspond to the tighter axes). 

While MOT2 is continuously on, MOT1 is operated in a pulsed regime. First, 
the trapping light and magnetic field are on for 100 ms and the MOT loads from the 
surrounding Rubidium vapour. Then, the light is switched off and a blue-detuned 
pushing laser beam, aligned on the vertical axis linking the two MOTs and passing 
through the differential pumping tube, is switched on for 6 ms. Because of the radiation 
pressure force, this light pulse transfers the atoms captured by MOT1 to MOT2 within 
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Pusher beam 




Figure 6. Schematic representation of the apparatus with the two magneto-optical 
trap chambers, the differential pumping tube, and the magnetic trap coils. The red 
arrows symbolize collimated laser beams. 

15 ms. The force exerted by the pushing beam is not strong enough to overcome 
the trapping force of MOT2. After typically 5 to 10 seconds of such loading, MOT2 
contains ~ 10 10 atoms and MOT1 is then switched off (light and magnetic field). The 
cloud is then compressed by a temporal dark MOT (compressed MOT): the cooling 
light detuning is changed from 5 = —3.5 r to 5 = —8 T and the magnetic field gradient 
is increased to B' 2 = 65.5 G/cm. This reduces the multiple-scattering-induced repulsive 
interaction between atoms and causes the cloud to shrink, thus increasing the density 
and collision rate by a factor of 3. The cloud is then further cooled to 80 /iK by a 
3-ms-long optical molasses phase (the field is switched off, and the detuning changed to 
5 = -10 r). 

For magnetic trapping, the atoms are then optically pumped to the |5 2 S'i/ 2 ,F = 
2, rriF = 2) Zeeman substate by a beam detuned by Szp = +3.2 T from the \5 2 Si/ 2 , F = 
2} — > \5 2 P 3 / 2 ,F = 2) transition and a repumper beam, detuned by 5zp rep . = — 3 T is 
applied to the \5 2 Si/2,F — 1) — > \5 2 Ps/2,F = 2) transition. A homogeneous magnetic 
field of ~ 0.5 G is aligned with the light wave vector to define the quantization axis. 
This optical pumping stage lasts 300 fis and then, all the light is switched off and a 
quadrupole magnetic field (54.1 G/cm) is abruptly turned on to trap the cloud. This is 
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followed by an adiabatic compression of the cloud, performed by linearly increasing the 
magnetic gradient to 278 G/cm within 500 ms. This compression stage is primordial 
to increase the elastic collision rate to a high enough value, an important requirement 
for evaporative cooling. At this point, the number of atoms is iV ~ 5 x 10 9 , and 
the temperature T ~ 400 /xK. In order to suppress Majorana losses, the quadrupole 
magnetic trap is then converted into a Ioffe-Pritchard trap by adiabatically ramping 
the current in a third coil (quadrupole-Ioffe-configuration or QUIC trap [23J) within 
500 ms. For cold enough atoms, this anisotropic trap is harmonic with radial and axial 
frequencies of 235.8 Hz and 22.2 Hz, respectively. Once the cloud is in the Ioffe-Pritchard 
trap, radio- frequency (rf) evaporative cooling is performed by ramping the rf frequency 
linearly from z/ start = 20 MHz to u stop = 1.3 MHz within 10 s. We are able to produce 
almost pure BECs (no discernible thermal fraction) containing up to 5 x 10 5 atoms. In 
order to produce an ultracold thermal cloud, the loading time of MOT2 is reduced to 
a few seconds. In this case, we are left, at the end of the evaporation, with a dilute, 
thermal gas, with a low elastic collision rate. 



3.1.2. Control of the trapping frequencies Implementing shortcuts to adiabaticity 
requires a precise control of the trapping frequencies, in a dynamical fashion. In our 
QUIC magnetic trap, this can be achieved by varying the current iq running through 
the 3 coils, and the current %b running through an additional pair of Helmoltz coils 
disposed along the long (axial) dimension of the trap (compensation coils). The resulting 
potential is 



U(x, y, z) = /u|B| ~ ji 



5o + ^(x 2 + , 2 ) + ^V 



(59) 



where fi/h = 1.4 MHz/G for atoms in \5 2 Si/2,F = 2 ) mp = 2). B' is the radial 
magnetic field gradient while B" corresponds to its curvature along y. The harmonic 



approximation of Eq. (59) describes accurately the potential seen by cold enough atoms 



i.e. ksT <C fiB [21]. Then, the radial and axial angular frequencies are 

H B'(i Q ) _ HI 



w ± =u z ~u s ~J , , U|| =u y =\—yJ B"(i Q ). (60) 

V m ^/B (i Q ,i Bo ) V m V 

These expressions show that the radial and axial frequencies can be controlled 
independently to some extent. The dynamical control of the currents is achieved using 
homemade, computer-controlled electronic circuits. The experimental realization of 
shortcut trajectories requires a careful preliminary calibration of the frequencies versus 
currents, which is achieved by monitoring the cloud's centre-of-mass oscillations after 
a small excitation. Due to the finite time response of the controlling circuit, it is also 
necessary to check the behavior of the frequency during an actual trajectory. This 
is illustrated in Fig. [7j where we compare the theoretical decompression trajectory of 
Fig. [2] (line) and measured experimental values (circles). In this example, the deviation 
is below 5%. 
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Figure 7. Vertical trap fre- 
quency calibration. The solid 
line is the theoretical shortcut 
decompression trajectory, the 
circles are the measured fre- 
quencies. The parameters of 
the decompression are given 
in Tab. El 



3.2. Shortcut to adiabaticity for a non-interacting gas 



For this first experiment, we use the procedure described in Sec. |3.1| to produce a thermal 
gas (N ~ 10 5 , T = 1.6 fiK) in the compressed trap of frequencies u x (0)/2n = 228.1 Hz, 
cu 3/ (0)/27r = 22.2 Hz and u z (0)/2tt = 235.8 Hz. The initial cloud-averaged collision rate 
per particle is 7 e i — 8 Hz, which corresponds to a collision time of ~ 125 ms. This is 
30 times the oscillation period, and more than 3 times the decompression time, which 
justifies the non-interacting approximation. 



We use here the decompression trajectory discussed in Sec. |2.1.2[ adapted to the 
vertical axis (Oz), with the parameters of Tab. [2j To maximize the decompression factor 
7 2 = lj z (0)/ 'cu z(tf), the compensation coils current ib is increased from iB (t = 0) ~ A 
to iB (tf) = 3.0 A, while the QUIC current is decreased from ig(t = 0) = 26.7 A to 
iq(tf) = 3.6 A (see the resulting trajectory in Fig. |7j. The decompression duration is 
tf = 35 ms. 

In theory, starting from a gas at equilibrium and temperature T in the compressed 
trap, a shortcut to adiabaticity should lead to an equilibrium state in the final trap, 
with a temperature Tf = T Q uj(tf)/u(0). This corresponds to a situation where entropy 
has not increased. On the contrary, for a non-optimal decompression, one expects to 
observe oscillations of the cloud's size and centre of mass in the decompressed trap, 
once the decompression is completed. To evaluate the efficiency of our shortcut, we 
thus perform the fast decompression, and hold the cloud in the decompressed trap for a 
variable amount of time. The trap is then abruptly switched off, and an absorption image 
is taken after a constant time of free expansion (6 ms). The amplitude of the dipole 
(oscillation of the centre of mass) and breathing modes (oscillation of the size) give access 
to the excess energy provided to the cloud, as compared to an adiabatic modification of 
the potential. If the cloud is reasonably at equilibrium after decompression, one can also 
directly measure the final temperature by measuring the evolution of the size during a 
free expansion. 

In the following, we compare four decompression trajectories: 
(i) the shortcut, given on Figs. [2fl and [f| 
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(ii) a linear decompression of the same duration (35 ms), 

(iii) an abrupt decompression, which, somehow, corresponds to a worst case scenario (in 
practice, the decompression time is 0.1 ms and u(t) is not controlled, it is imposed 
by the response of the magnetic trap control electronics), 

(iv) a 6-s-long linear decompression, which can be considered nearly adiabatic. 

What is referred to as 'linear decompression' corresponds to both control currents being 
varied linearly with time. The corresponding frequency trajectory is not linear. 

The experimental results are summarized on Fig. [8] In the case of the 6-s-long 
linear ramp (filled squares), very little residual excitation is observed (although the 
residual dipole mode is still measurable). In the shortcut case (open circles), clear 
oscillations of the cloud width and centre-of-mass position are seen, but they are much 
reduced compared to the fast linear ramp (diamonds) and abrupt decompression (open 
squares) . 




Figure 8. Comparison between different trap decompression schemes (along the 
vertical axis). Open red circles: shortcut decompression in 35 ms; black diamonds: 
linear decompression in 35 ms; solid blue squares: linear decompression in 6 s; open 
black squares: abrupt decompression. The decompression is performed, and then, the 
cloud is held in the decompressed trap for a variable time. We monitor (a) the vertical 
centre-of-mass position (dipole mode) and (b) the cloud size (breathing mode), after 
6 ms time of flight. On (a), the solid lines are sine fits, on (b) they just connect the 
points to guide the eye. 

Compared to the linear decompression in 35 ms, the shortcut reduces the amplitude 
of the dipole mode by a factor of 7.2 (obtained from the sine fits) and the amplitude 
of the breathing mode by a factor of 3 (comparison of the standard deviations of the 
two sets of data). The excess energy, which is dominated by the centre-of-mass energy, 
is thus reduced by a factor of ~ 52. In the case of the 6-s-long ramp, we measured a 
final temperature of the cloud of 130 nK, a factor 12.5 below the initial one. This is 
consistent with the expected value of 15. The small difference may arise from a small 
heating rate due to the fluctuations of the magnetic trap. 
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The fact that the shortcut decompression still produces sizeable excitations is due 
to experimental imperfections. Several possible causes can be invoked. Firstly, as seen 
on Fig. [7j there are still small deviations from the ideal trajectory. These may have an 
impact, especially in the last phase of the trajectory where the cloud is subject to a large 
acceleration (see Fig. [3]). Second, as can again be seen in Fig. |3j during the trajectory 
the cloud wanders quite far (several hundred /xm) from the trap centre and feels the non- 
harmonic part of the potential. This effect is difficult to quantify since our knowledge 
of the potential shape is not sufficiently accurate (however, the anharmonicity could be 
inferred from variations of the oscillation frequency with amplitude). In principle, it 
could be avoided by designing other shortcut trajectories keeping the cloud closer to the 
trap centre at all times. 




40 80 120 160 200 
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Figure 9. Impact of the ver- 
tical decompression schemes 
on the axial size (y direction) . 
Same colors and symbols as in 
Fig. [8] The amplitude of the 
axial breathing mode is not 
affected by the use of a short- 
cut trajectory adapted to the 
radial dimensions. 



Fig. [9] shows the behavior of the axial size of the cloud in the conditions of Fig. [8Jd. 
Since the shortcut trajectory was designed only for the radial dimensions, the resulting 
axial breathing mode is of the same magnitude as for the linear decompression. 

We compare on Fig. 10 the results of the shortcut decompression to linear ramps of 
various durations. The vertical axis in this figure represents amplitudes of oscillations 
after trap decompression, either of the centre-of-mass position (filled symbols) or of 
the cloud radius (open symbols), scaled by their values for an abrupt decompression 
(tf ~ 0.1 ms). The horizontal axis is the duration of the decompression tf (notice the 
logarithmic scale). The circles correspond to linear decompressions while the stars are 
the shortcut results. As can be seen, fulfilling the adiabaticity criterion is easier for the 
breathing mode (size oscillation) than for the dipole mode (centre-of-mass oscillation): 
the oscillation amplitude is reduced by a factor of 2 for tf = 20 ms for the earlier, and 
for tf ~ 150 ms for the latter. Using the amplitude of the dipole mode as a criterion 
to compare the linear and shortcut schemes, one sees that the decompression time is 
reduced by a factor of 37. 
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Figure 10. Comparison between linear and shortcut decompression schemes. We 
plot the scaled oscillation amplitudes of the breathing (cloud size, open symbols) and 
dipole (centre-of-mass position) modes versus the decompression time tf. The circles 
and stars correspond to linear and shortcut decompressions, respectively. 



3. 3. Shortcut to adiabaticity for an interacting condensate 

As opposed to the previous case of non-interacting atoms, the decompression of a BEC 
is an intrinsically 3D problem because of the interactions. As a result, both the radial 



and axial frequencies have to be varied following Eqs. (56) and (57) in order to realize a 
shortcut to adiabaticity. In the present section, we describe a decompression experiment 
based on the trajectories discussed in Sec. 2.2.1 and represented in Fig. |4j In this scheme, 



the radial frequency is decreased by a factor of 9, while the axial frequency is adjusted 
to maintain the axial size of the BEC fixed during the whole trajectory. Accordingly, 
the axial frequency is decreased by a factor of 3. 

We start from an initial BEC containing 1.3 x 10 5 atoms in the condensed fraction, 
and 7 x 10 4 non-condensed atoms at a temperature of 130 nK. The experimental scheme 
is similar to that employed for the thermal cloud. Here, we use a longer time of flight of 
28 ms to characterize the various excitations generated by rapid decompressions. Three 
decompression schemes are compared: 

(i) the shortcut to adiabaticity in 30 ms, 

(ii) the linear decompression in 30 ms, 

(iii) an abrupt decompression. 

Contrary to the previous case of a thermal cloud, the BEC cannot be held for more 
than 150 ms in the compressed magnetic trap because of a relatively high heating rate. 
Thus, here we cannot compare our scheme to the adiabatic limit corresponding to a 
slow linear decompression. 



Figure 11 shows the temporal behaviour of the cloud following the linear and 
shortcut decompressions. These absorption images are taken in the (y, z) plane, after a 
certain holding time in the decompressed trap (indicated in the figure) plus a 28-ms-long 
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time of flight. The field of view is 545 /iin x 545 /im. The centre-of-mass motion has been 
subtracted from these data for better clarity. In the linear case the BEC (yellow central 
part) experiences large deformations and oscillations of its aspect ratio, whereas in the 
shortcut case it remains nearly perfectly stationary. Surprisingly, in the case of the linear 
decompression the BEC also oscillates angularly. We attribute this to an uncontrolled 
tilt of the trap axes during the decompression. This will be discussed in more details 
later. The nearly isotropic aspect of the BEC after the shortcut decompression is due 
to the value of the time of flight, which is close to the critical time of inversion of the 
aspect ratio. The thermal component surrounding the BEC (red halo) is also visible. 
It's temporal evolution is discussed at the end of this section. 



a) 




• 


• 


■0fr r 


40^' 








20 ms 


30 ms 


40 ms 


50 ms 


60 ms 


70 ms 


80 ms 


90 ms 




















isity 


ft 








mi- 






I 


der 


100 ms 


110 ms 


120 ms 


130 ms 


140 ms 


150 ms 


160 ms 


170 ms 


ical 




| Opt! 


b) 
















• 


• 


• 




m 






m 




20 ms 


30 ms 


40 ms 


50 ms 


60 ms 


70 ms 


80 ms 


90 ms 




# 


• 




m 


m 


• 




• 




100 ms 


110 ms 


120 ms 


130 ms 


140 ms 


150 ms 


160 ms 


170 ms 





3 



0.1 



Figure 11. Comparison of linear and shortcut BEC decompressions. We compare the 
time evolution of the BEC after two different decompression schemes: (a) a 30-ms-long 
linear ramp and (b) the shortcut trajectory (see text). The centre-of-mass motion has 
been subtracted from these time-of-flight images for clarity. On each image, the region 
where the optical density is highest (yellow and white) correspond to the condensate, 
while the red halo is the thermal component. 



To provide a more quantitative analysis, the column densities obtained from the 
absorption images were fitted with a 2D bimodal distribution consisting of a Gaussian 
component, accounting for the thermal fraction, plus a 3D inverted parabola integrated 
along one dimension, accounting for the condensed atoms. The fitting parameters were 
the cloud centre, two angles, one for each couple of eigenaxes of each components, and 
the two widths of each components. 

In Fig. 12 i) is reported the centre-of-mass oscillations (dipole mode) for the abrupt 
(squares), linear (diamonds) and shortcut (circles). Figure [12)3) shows the oscillations of 
the BEC's aspect ratio (breathing mode). All measurements are performed after a 28 ms 
time of flight. As in the case of the non-interacting cloud, the shortcut scheme reduces 
the amplitude of the dipole mode compared to a standard linear decompression, here by 
a factor of 4.3. For our relatively long time of flight, the measured positions reflect the 
atomic velocities. Thus, using the shortcut scheme reduces the kinetic energy associated 
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with the dipole mode by a factor of 18.5 compared to the linear one (and 36 compared to 
the abrupt). The residual energy after the shortcut decompression is 580 nK. As can be 



seen in Fig. 12 d), both non-optimal schemes induce very large oscillations of the BEC's 
aspect ratio, with a rather complicated dynamics. A Fourier analysis reveals a main 
oscillation frequency of 47 Hz, consistent with a radial breathing mode at 2cu_i_ [2"oT42"7j . A 
smaller contribution at 12.5 Hz corresponds to the axial breathing mode at \Jh/2uj\\ |27j . 
The shortcut scheme suppresses strikingly these breathing oscillations, yielding a BEC 
close to the targeted equilibrium state. 
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Figure 12. Decompression-induced excitations of the BEC. We report the temporal 
evolution of (a) the centre-of-mass position and (b) the aspect ratio of the BEC after 
three different decompression schemes: an abrupt decompression (black squares); a 
30 ms linear ramp (black diamonds); the 30 ms shortcut trajectory (red circles). All 
measurement are performed after 28 ms of time of flight. 



As emphasized in section [272] , the shortcut trajectory employed in this experiment is 
also valid for the thermal fraction, in the radial dimensions only. This is demonstrated 



in Fig. 13, where we compare the oscillations of the radial (open symbols) and axial 
(filled symbols) sizes of a) the BEC, and b) the thermal fraction, after the shortcut 
decompression. The BEC's TF radius is stationary with an average value of 46.8 fivcv 
close to the theoretical value (43/xm). As can be observed in Fig. fl3b), the radial size 
of the thermal fraction is also quite stationary as expected from a shortcut trajectory. 
Thus, this experiment demonstrates that both a non-interacting thermal gas and an 
interacting BEC can be decompressed simultaneously using an appropriate shortcut 
trajectory. The observed behavior is also qualitatively consistent with our initial 
assumption that the BEC and thermal fraction are independent. However, we expect 
that ultimately the validity of this approach will be limited by the interaction between 
the condensed and non-condensed fractions. The temperature inferred from the radial 
size of the thermal component is 22 nK, a factor of 6 below the initial one. This 
factor is smaller than the expected one (oo±(0)/uj±(tf) = 9), and even if we improve 
the experimental set-up to realize the ideal frequency trajectory we would probably be 
limited by the transfer of energy from the axial breathing mode via the interaction with 
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the condensate. Indeed, the axial size of the thermal fraction presents clear breathing 
oscillations, reflecting the fact that the shortcut trajectory u» (t) is not valid in this case, 
as expected. 
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Figure 13. BEC versus thermal cloud decompression. We plot (a) the sizes of the 
BEC and (b) thermal component versus the time spend in the decompressed trap for 
the shortcut trajectory. The filled and empty symbols correspond to the axial and 
radial (vertical) directions respectively. The line is a sine fit to the thermal fraction 
axial size. 



A striking feature in Fig. [TT^ i was the large angular oscillation of the BEC after 
the linear decompression. This unexpected effect is due to a slight tilt of the QUIC 
trap eigenaxes (3°) in the (y, z) plane as the trap centre moves downwards due to 
gravity. Because of this, an angular momentum is imparted to the atoms during the 
decompression, exciting a 'scissors mode' [28] [29]. Our nearly critical time of flight then 



results in a magnification and a deformation of the scissors oscillations [301 131] . Fig. 14 
shows an example of these oscillations, together with a GPE simulation (red line). 



4. Other possible applications 

In this section, we attempt to generalize the shortcut decompression of Bose-Einstein 
condensates to other situations which may find applications in experiments where a 
fast and large modification of the width of the velocity distribution or of the chemical 
potential is required. 



4-1. Arbitrary variation of a harmonic potential 

Let us consider the time evolution of a condensate in the time-dependent harmonic 
potential of the form 

U(r, t) = r*W(t)r + r*u(t) (61) 
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Figure 14. Experimental observation of a scissors mode excitation following the linear 
decompression (diamonds). The red line is a GPE simulation. The oscillation is not 
quantitatively reproduced because it depends on the precise way the trap is rotated 
during decompression, which is not known for the whole trajectory. Only the final 
tilt of 3° was measured. For the GPE simulation, the trap angle was assumed to be 
proportional to the trap bottom displacement from its original position. 



where the symmetric matrix W(t) = R 1 (t)W(t)R(t) represents the harmonic potential 
of stiffness 

/ u*(t) \ 
W= wj(t) , (62) 
V u 2 z (t) J 

rotated by a rotation matrix R(t). The column vectors r and u respectively represent the 
position and a spatially homogeneous force which may depend on time. The superscript 
* indicates the transpose of vectors or matrices. 

To solve Eq. (31) we look for a linear change of variables p(r, {^(t)}, {rf m (t)}) 
where the bi/s are scaling and rotation functions for the r^'s. Let B be a 3x3 matrix 
which elements are the functions by. The transformation is 



p = B-\t) (r - r cm (t)) = £ _1 (t)r + a(t). 



(63) 

In the TF limit, and if the matrix BB~ l is symmetric, Eq. ( |31~j ) is invariant under this 
transformation. The full derivation is given in Appendix A, but we give here the key 
elements. 

The TF approximation consists in neglecting the kinetic-energy-like term 

d 2 x 



i,j,k 



dpidp k '' 



(64) 



x(p,t) being defined as in Eq. (32). In this regime, the condensate wavefunction 
x(p, T ) verifies the equation of motion Eq. (38 ), under the action of the time-independent 
potential 



U(p,0) = -mp*iy(0)p 



(65) 
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if the generic scaling functions satisfy 



B l B + B l WB = V P^-, (66) 
det-D 

r cm + W(t)r cm - — u = 0. (67) 
m 

It is worthwhile recalling that, as shown by the above equations, the evolution of B 
is decoupled from the centre-of-mass motion which evolves with the net external force. 
The phase of the wavefunction is chosen as 

m \ 1 1 

«Hr,f) = j < -r'BB-'r - r'Bd\ + «l„(f), (68) 

jU<(a'B>Ba-a>-^a). (69) 



with 

The wavefunction normalization is 



and the time r is defined by 



A = (detfi)- 1 / 2 , (70) 
dT 1 (71) 



dt detB 



The derivation of the scaling equations (Appendix A) relies on the particular choice of 
the above phase <ft which verifies 

V r = -^5^ or v(r) = BB- 1 r-BB- 1 r cm + B- 1 t cm , (72) 
n dt 

v(r) being the velocity field of the condensate, and on the assumption that the matrix 
BB^ 1 is symmetric. The first condition consists in imposing that there are no terms 
linear in momentum in the GPE in the p-coordinate frame; if the first condition is 
fulfilled the second imposes that the velocity field is irrotational, namely that the 
condensate is a superfluid everywhere as well. This implies that our scaling ansatz 
does not take into account the presence of quantized vortices and thus can describe the 
dynamics of a rotated condensate only below the critical angular velocity a c ~ O^cu^ for 
a slightly anisotropic confinement [32J, or in general, for a metastable configuration [33J. 
Nevertheless, a slightly modified ansatz could be deviced to incorporate the possibility 
of quantized vortices. It is also possible to relax the first condition and allow for 
terms in the GPE that contain for instance the angular momentum components. These 
extensions are deferred for future studies. 



Equations (66) and (67) can be used to determine the dipolar, compressional 



and scissors modes for a harmonically-trapped superfluid condensate (see Appendix 



[b| . Replacing detS with (detB) 13 in Eq. (66), the same equation describes the 
compression and the scissors dynamics of a superfluid characterized by an equation 
of state n(n) oc n", as it has been already shown for the quadrupolar modes [31] and as 
it can be easily deduced by Eq. (A. 8) of the Appendix. In the following we present three 
possible shortcut trajectories based on these scaling equations and adapted to compress 
or decompress and rotate a BEC in the absence and in the presence of gravity. 
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4-2. Uniform decompression or compression of a condensate 

We now consider the particular case of u = and W diagonal. If one wants 
to compress or decompress the condensate without modifying the condensate aspect 
ratio, the condition Ui(tf) = Co>j(0)/7 2 must hold for any i. The boundary conditions 
for the shortcut solution are: bu(0) = bu(tf) = 0, 6^(0) = 1, &«(£/) = 7 4//5 and 
bu(0) = buitf) = 0. One possible solution is to set all 6ji(t)'s equal to a unique function 



k=0 



with c = 1, ci = c 2 = 0, c 3 = 10(7 4 / 5 - 1), c 4 = -15(7 4 / 5 - 1), c 5 = 6(7 4 / 5 - 1). The 
time evolution of the trap frequencies C0i(t) will be given by the equation 

If the kinetic energy is negligible during the whole decompression, the final state is a 
BEC at equilibrium with a chemical potential that has been divided by a factor of 7 16 / 5 
(because ji oc (IljWj) 2 ^ 5 ). 

4-3. General compression or decompression in the presence of gravity 

We now consider the case where W(t) is diagonal with u x (t) = u z {t) = u)±(t), 
u) y (t) = w\\(t), and u z = rag. A general compression or decompression of a condensate 



confined in this axially-symmetric trap (51) can be realized in two steps: (i) in the first 



step {t G [0, £]), b\\ is kept fixed as outlined in Sec. 2.2, while the desired final value of 
b± = b±(tf) (R±(tf)) is reached; (ii) then (t £ [i,tf]) b± is fixed and b\\ evolves according 
to the set of equations: 

*® ^ W' (78) 

h(t) + b ll (t)uf l (t) = ^, (76) 
h(t)c(t) = ul(t) (c(t) - b {] (t)) , (77) 



where c(t) = —uj_(t)r c z m (t) / (gb±(t)) as in Eq. (52). Also in this case one can write the 
function c(t) as a polynomial of order > 9 (see Eq. (50)) with the first coefficient fixed 
to one and the following four coefficients fixed to zero. The other coefficients are fixed 
by the boundary conditions at the time tf of the function c(t) and of the function b»(t), 



that from Eq. (77) can be written as 



ul(t)c(t) 

6||(t) -"c(t)-a;i(t)' (7f 
and by the boundary conditions of their derivatives at the same time tf. 
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4-4- Rotation of the BEC in the presence of gravity 



Now we propose a shortcut trajectory to rotate an axially-symmetric BEC of an angle 
a, in the presence of the gravity In this case 



/ wi(0) 



W(0) = 

and W(t f ) = R^W{0)Ra, with 

Ra z 







wj?(0) 



\ 


.(0) ) 



(79) 



/ 1 



cos a sin a 
— sin a cos a 



(80) 



Let us suppose, for instance, w_l(0) < W||(0), with a;ii(0) = Xu±(0). The tilted ground- 
state for the potential W(tf) can be obtained in two steps: (i) during a time i, fixing 
fey, decompressing the BEC in the radial direction up to the value b±(t) = A -1 . At t = i 
the trap is spherical with frequency Q = Xu\\(0) and the BEC is spherical with a TF 
radius equals to -R||(0). (ii) Fixing b» along the direction y', compressing in the direction 
x' and z', where the axis r' are defined by r' = i? a r. Using the new coordinate reference 
frame, and setting c z (t) = — u 2 rl m (t) / (gb±(t) cosa), and c y {t) = — a) 2 r^ m (t)/(5fsina), 
we obtain the set of equations 



b ± (t) + b ± (t)ul(t)=Q 2 /b 3 ± (t), 
u\\(t) = ti/b ± (t), 
b\{t)c z {t) + 2bl(t)b ± (t)c z (t) + u 2 (c z (t) - bl(t)) = 0, 
bl(t)c y (t)+Cj 2 (c y (t)-b\{t)) =0, 



(81) 
(82) 
(83) 
(84) 



the latter describing the centre-of-mass motion in the y' direction. The boundary 
conditions for such a problem are: b±(i) = c z (t) = c y (t) = 1, b±(tf) = A, c z (tf) = A 3 , 
c y (tf) = A 2 , and that all the first and the second derivatives with respect to time are 
null at t — i and tf. In this case a finite-order polynomial ansatz in r for q was found 
to be inadequate as a solution of the scaling equations due to the coupling of c y and 
c z . A full numerical solution of the dynamical equation using, e.g., a shooting method 
[35] or following a strategy as that implemented in optimal control [9] may be needed 
in finding a shortcut trajectory in this case. 



5. Conclusion 



We have experimentally demonstrated the controlled transfer of trapped ultracold atoms 
between two stationary states using a faster-than-adiabatic process which reduces the 
transfer time down to a few tens of milliseconds. The transfer is achieved by engineering 
specific trajectories of the external trapping frequencies that take explicitly into account 



Shortcuts to adiabaticity for trapped ultracold gases 



30 



the spatial shift introduced by gravity. This scheme was successfully applied both to 
a thermal gas of atoms and to an almost pure Bose-Einstein condensate. The scheme 
used is flexible enough to be adapted to both situations even though in the thermal 
gas interactions does not play a significant role while the Bose-Einstein condensate is 
strongly affected by the s-wave scattering of atoms. The residual excitations observed 
after the shortcut decompressions in the present demonstration experiments are due to 
our imperfect control over the time-varying magnetic trapping potential, and could be 
substantially reduced in future realizations. 

Theoretically, the design of the transfer process was based on the invariants of 
motion and scaling equations techniques which turned out to be possible thanks to 
the harmonic shape of the external potential. In our scheme, the invariant of motion 
technique (for non-interacting particles) and the scaling equations technique (valid for 
both the non-interacting and the interacting gas) are tightly connected. The invariant 
of motion we used is a time-independent harmonic oscillator Hamiltonian that can be 
obtained by a time-dependent canonical transformation of position and momentum. In 
the scaling equation technique, we looked for a scaling plus shift transformation of the 
coordinate that allowed the equation of motion for the system to be time-independent 
(except for terms that are not coordinate-dependent). In both cases the whole dynamics 
is included in the new set of (canonical) coordinates, that depend on the trap frequencies. 
We also showed that these techniques can be generalized to include the rotation of the 
eigenaxes without much effort. 

Very often, in cold-atom experiments, samples are prepared by transferring atoms 
from some confinement to another, e.g., from a magneto-optical trap to a magnetic 
quadrupolar trap, from a quadrupolar trap to a Ioffe-Pritchard trap, from an harmonic 
confinement to an optical lattice, etc., the main limitation being that, for short transfer 
times, parasitic excitations may show up. The main application of our scheme is to 
guide this transfer in order to prepare a very cold sample in a very short time with the 
desired geometry and without exciting unwanted modes. The shortcut-to-adiabacity 
scheme proposed here could be applied to non-interacting particles such as cold gases 
or ultracold spin-polarized fermions, to normal or superfluid (bosonic or fermionic as 
well) gases in the hydrodynamic regimes, and to strongly correlated systems such as 
the Tonks gas. In this paper we focused on explicit solutions to transfer atoms between 
two stable states, but the same strategy could be applied to control the generation of 
metastable states, vortex states, or some exotic out-of-equilibrium states. We plan to 
explore these possibilities in future studies. 

Acknowledgments 

This work was supported by CNRS and Universite de Nice-Sophia Antipolis. We also 
acknowledge financial support from Region PACA, Federation Wolfgang Doeblin, and 
CNRS-CONICET international cooperation grant n°22966. JFS acknowledges support 
from the French ministry of research and education for his funding, and thanks Mario 



Shortcuts to adiabaticity for trapped ultracold gases 



31 



Gattobigio and Michel Le Bellac for helpful discussions on theoretical aspects. 



Appendix A. Demonstation of Eqs. (66) and (67) 



In this appendix we derive Eqs. (|66|)-(|71|). The starting point is the GPE (31) for a 

(A.l) 



general potential (61). We look for a solution of the form 

^(r,t)=A(t) X (p,T)e^ 
with 

P 



(A.2) 



Equation (31) then takes the form 



ih 



A _ dp(B,a) dxdr ■ 



2m 



kj 



d 2 X 
dpidp k 



+ 2i(B~ 1 V r 4>) ■ V pX + i{Vl<l>)x ~ (V r ) 2 x \ (A ' 3) 



+-m { [B(p - a)fW[B{p - a)}} X + u t B{p - a) X + g\A\ 2 \ X \ 2 X. 



We look for the conditions that A, B, and a have to verify aiming to simplify Eq. (A.3) 
to the form 



ih^xiPiT) = [u(p,0) + VN\x(p,r)\ 2 ] X (p,r), 



(A.4) 



in the TF limit, namely, neglecting the kinetic term given in Eq. (64). We deduce 
immediately that (i) the second term of Eq. (A.3) has to be equal to the sixth, and (ii) 
the first to the seventh. Condition (i) leads to 

_ m 



h 



BUB-^r + a 



(A.5) 



that has a solution if the matrix B[B x ] = —BB 1 is symmetricjlj If this condition 
holds, we get Eq. (68) for the phase 0. Condition (ii) can be written as 

1 



AA 



-i 



-tr(BB 



(A.6) 



Using the invariance of the trace and of the determinant, the evolution of A can be 
rewritten in term of the eigenvalues /V of the matrix B as 



AA 



-i 



d 



dt 



\nA 



Id 
2dt 

ld_ 
2dt 



IndetS 
In detS 



(A.7) 



| In a general case the matrix BB -1 can be split into a symmetric and an antisymmetric part. In the 
p-frame of reference, the antisymmetric part gives rise to a rotational term proportional to the angular 
momentum and only the symmetric part of BB^ 1 contributes to the phase of the wave function. The 
rotational term can be neglected for nearly-isotropic trap or for small angular velocities of the trap. 
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If, e.g., at t = we have that B is the identity and A — 1, equation (A. 7) yields Eq. (70). 

Moreover from the comparison between the third term in Eq. ( |A.3 ) and the non- 
linear term (condition (hi)), we deduce Eq. (71). Taking into account (i)-(iii), Eq. (A. 3) 
reduces to 



ih 



dx 



dr 
detB 



MetB 



m 
~2 



dt 

BB~ l v 



Bit 2 + ^BB-h + r'BiB-^r 



-m^Ba - mr*Bd + -m { [B{p - a)]*W[£(p - a)] } + w*5(p - a) 

+ g\x\ 2 x- 



(A.8) 



X 



By imposing the quadratic term in p to be equal to yp'iy p, we get condition (iv), 
i.e., Eq. (66); the fifth condition is that the linear term in p vanishes and thus leads to 
(67); finally by requiring that the p-independent term be null, we get (69) for 4> . 



Appendix B. Low-lying modes 



Equation (67) describes the dipolar modes for the centre of mass and Eq. (66) the 
quadrupolar and the scissors modes. The low-lying eigenfrequencies of these latter 
modes can be obtained by solving the equation of motion for the matrix B for the case 
of a tilt of the trap of a small angle a. At t > 0, the matrix W is constant and can be 
written as 

/ oo 2 



W 








u f\ a ( u f\ — w i) 



W° + 6W 



where 



a(u 2 - wl) 
( 



and 



5W 



( 




V 



< 











"? 











CO 2 








V 



a(u 2 — bj\ 








/ 



We look for solutions of the form B l = 1 + 5. Equation (66) takes the form: 

5 ~ -W°5 - 5 l W° - (Tt5)W° + 5W, 
at the first order in 5. For the diagonal terms we have 

Su = -2u 2 5u - (Tr5)u 2 . 



(B.1) 



(B-2) 



(B-3) 



(B-4) 



(B.5) 
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Setting 5a = Aie lQ,t , we obtain the following coupled equations 

-Q 2 A X = -2u 2 ± A x - (A x + A y + A 2 )wi 
-fi 2 A y = -2u|A y - (A x + A y + A z )u\ (B.6) 
-H 2 A Z = -2ulA z - (A, + A y + A 2 )ur[ 

whose solutions are the surface mode VL = y/2u± for any values of u± and and the 
breathing modes f2 ~ 2co± and ~ y/h/2oj\\ in the cigar-shape regime <C wi. 



For the off-diagonal terms 5^, = {2,3} or {3,2}), Eq. (B.4) gives 



Sij = -u 2 Sij - UjSji + a(u 2 - u 2 ± ) (B.7) 

namely 

Sij + 5ji = -(oo 2 + wj)(<% + <y + 2a(u\ - wi). (B.8) 



Equation (B.8) has solution 



5 23 = $32 = a 11 n2 X [1 - cos(fit)], (B.9) 

where £l s = (lo^+ui 2 ) 1 ^ 2 . This is a scissors mode with boundary conditions Sij(t = 0) = 
and 6ij(t = 0) = 0. 
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